#drafted by Jeremy VanDerWal ( jjvanderwal@gmail.com ... www.jjvanderwal.com )
#GNU General Public License .. feel free to use / distribute ... no warranties

################################################################################
#load the library
library(SDMTools)

#define if this is to be run @ UQ
UQ=FALSE

#define the working directory
base.dir = ifelse (UQ,'/home/uqvdwj/','/home/jc148322/')#'/ctbccr/scratch/') 
base.data.dir = paste(base.dir,'Climate/AWAP.direct.download/',sep='')
work.dir = paste(base.data.dir,'summaries/Oz/monthly.asc/',sep=''); dir.create(work.dir,recursive=TRUE); setwd(work.dir)

###########################################################################################
###########################################################################################
#create the base ascii grid file and the row/col of the pos data

#create the base asc
tasc = read.asc.gz(paste(base.data.dir,'rain/monthly/201101.asc.gz',sep=''))
#readin the base positions
pos = read.csv(paste('/home/jc165798/','Climate/OZmask/mask.pnts.csv',sep=''))
pos$order = 1:nrow(pos)
tpos = as.data.frame(which(is.finite(tasc),arr.ind=TRUE))
tpos$lat = getXYcoords(tasc)$y[tpos$col]; tpos$lon = getXYcoords(tasc)$x[tpos$row] #convert to lat & long
pos = merge(pos,tpos) #add row col information
pos = pos[order(pos$order),]
pos$order = NULL
#change all values of tasc to NA and then only fill in tasc values at positions to 1
tasc[,] = NA
tasc[cbind(pos$row,pos$col)] = 1
write.asc.gz(tasc,paste(base.data.dir,'summaries/Oz/base.asc',sep='')) #write out the asc file
write.csv(pos,paste(base.data.dir,'summaries/Oz/base.positions.csv',sep=''),row.names=FALSE)

###########################################################################################
###########################################################################################
#start creating monthly summary information

#define years/months of interest and variables of interest
yois = 2011:2012
mois = sprintf('%02i',1:12)
vois = c('bc01','bc04','bc05','bc06','bc12','bc15','bc16','bc17')
tois = c(2,5,11,35,59,119) #representing the preceeding 3,6,12 months

#read in the base asc & positions
base.asc = paste(base.data.dir,'summaries/Oz/base.asc.gz',sep='')
pos = paste(base.data.dir,'summaries/Oz/base.positions.csv',sep='')

#go to the temporary pbs directory
pbs.dir = paste(base.dir,'tmp.pbs/',sep=''); dir.create(pbs.dir); setwd(pbs.dir)

#create the output directories
for (yy in yois) { for (mm in mois) { dir.create(paste(work.dir,yy,mm,sep='')) } }

#cycle through the times
for (ttime in tois) {
	#cycle through the variables
	for (voi in vois) {
		if (!(any(c('bc16','bc17') == voi) & ttime < 11)) { #do not do anything less than annual for bc16 & 17
			##create the sh file
			zz = file(paste(voi,ttime,'.sh',sep=''),'w')
				cat('##################################\n',file=zz)
				cat('#!/bin/bash\n',file=zz)
				cat('source /etc/profile \n',file=zz)
				cat('module load R \n',file=zz) #load the necessary module				
				cat('cd $PBS_O_WORKDIR\n',file=zz)
				cat('R CMD BATCH --no-save ',voi,ttime,'.R\n',sep='',file=zz)
				cat('##################################\n',file=zz)
			close(zz)
			#create the R summarizaing script
			zz = file(paste(voi,ttime,'.R',sep=''),'w')
				cat('library(SDMTools)\n',file=zz)
				cat('#define directories\n',sep='',file=zz)
				cat('work.dir = "',work.dir,'"; setwd(work.dir)\n',sep='',file=zz)
				cat('#read in necessary data\n',file=zz)
				cat('base.asc = read.asc.gz("',base.asc,'")\n',sep='',file=zz)
				cat('pos = read.csv("',pos,'")\n',sep='',file=zz)
				cat('#define years/months of interest\n',file=zz)
				cat('yois = 1950:2010; mois = sprintf("%02i",1:12)\n',file=zz)
				cat('#bring in the primary datafile\n',file=zz)
				if(voi=='bc01' | voi=='bc04') {
					cat('tdata = read.csv("',base.data.dir,'summaries/Oz/monthly.csv/tmean19402010.csv")\n',sep='',file=zz)
				} else if (voi=='bc05') {
					cat('tdata = read.csv("',base.data.dir,'summaries/Oz/monthly.csv/tmax19402010.csv")\n',sep='',file=zz)
				} else if (voi=='bc06') {
					cat('tdata = read.csv("',base.data.dir,'summaries/Oz/monthly.csv/tmin19402010.csv")\n',sep='',file=zz)
				} else {
					cat('tdata = read.csv("',base.data.dir,'summaries/Oz/monthly.csv/rain19402010.csv")\n',sep='',file=zz)
				}
				cat('#ensure cols are in order of date\n',file=zz)
				cat('tdata = as.matrix(tdata); tdata = tdata[,order(colnames(tdata))]\n',file=zz)
				cat('#define some functions\n',file=zz)
				cat('new.sd = function(x) { tsum = sd(x,na.rm=T); if (is.finite(tsum)){ return(tsum) } else { return(0) } }\n',file=zz)
				cat('new.max = function(x){ return(max(x,na.rm=T)) }\n',file=zz)
				cat('new.min = function(x){ return(min(x,na.rm=T)) }\n',file=zz)
				cat('#cycle through each of the years / months of interest\n',file=zz)
				cat('for (yy in yois) {\n',sep='',file=zz)
				cat('    cat(yy)\n',sep='',file=zz)
				cat('    for (mm in mois) {\n',sep='',file=zz)
				cat('        cat(".")\n',sep='',file=zz)
				cat('        coi = paste("Ym",yy,mm,sep=""); coi.in = grep(coi,colnames(tdata)) #define the column of interest\n',sep='',file=zz)
				if(voi=='bc01') {
					if (ttime<12) { 
						cat('        tt = round(rowMeans(tdata[,coi.in:(coi.in-',ttime,')],na.rm=T),1) #extract the summary data\n',sep='',file=zz)
					} else { 
						cat('        tt = tdata[,coi.in:(coi.in-',ttime,')] #extract a subset of the data\n',sep='',file=zz)
						cat('        for (ii in 1:12) {tt[,ii] = rowMeans(tt[,seq(ii,ncol(tt),12)],na.rm=TRUE)}; tt = tt[,1:12] #summarize to get averages\n',sep='',file=zz)
						cat('        tt = round(rowMeans(tt,na.rm=T),1) #extract the summary data\n',sep='',file=zz)
					}
				} else if (voi=='bc04') {
					if (ttime<12) {
						cat('        tt = round((apply(tdata[,coi.in:(coi.in-',ttime,')],1,new.sd))*100,0) #extract the summary data\n',sep='',file=zz)
					} else { 
						cat('        tt = tdata[,coi.in:(coi.in-',ttime,')] #extract a subset of the data\n',sep='',file=zz)
						cat('        for (ii in 1:12) {tt[,ii] = rowMeans(tt[,seq(ii,ncol(tt),12)],na.rm=TRUE)}; tt = tt[,1:12] #summarize to get averages\n',sep='',file=zz)
						cat('        tt = round((apply(tt,1,new.sd))*100,0) #extract the summary data\n',sep='',file=zz)
					}
				} else if (voi=='bc05') {
					if (ttime<12) {
						cat('        tt = apply(tdata[,coi.in:(coi.in-',ttime,')],1,new.max) #extract the summary data\n',sep='',file=zz)
					} else { 
						cat('        tt = tdata[,coi.in:(coi.in-',ttime,')] #extract a subset of the data\n',sep='',file=zz)
						cat('        for (ii in 1:12) {tt[,ii] = rowMeans(tt[,seq(ii,ncol(tt),12)],na.rm=TRUE)}; tt = tt[,1:12] #summarize to get averages\n',sep='',file=zz)
						cat('        tt = apply(tt,1,new.max) #extract the summary data\n',sep='',file=zz)
					}
				} else if (voi=='bc06') {
					if (ttime<12) {
						cat('        tt = apply(tdata[,coi.in:(coi.in-',ttime,')],1,new.min) #extract the summary data\n',sep='',file=zz)
					} else { 
						cat('        tt = tdata[,coi.in:(coi.in-',ttime,')] #extract a subset of the data\n',sep='',file=zz)
						cat('        for (ii in 1:12) {tt[,ii] = rowMeans(tt[,seq(ii,ncol(tt),12)],na.rm=TRUE)}; tt = tt[,1:12] #summarize to get averages\n',sep='',file=zz)
						cat('        tt = apply(tt,1,new.min) #extract the summary data\n',sep='',file=zz)
					}
				} else if (voi=='bc12') {
					if (ttime<12) {
						cat('        tt = round(rowSums(tdata[,coi.in:(coi.in-',ttime,')],na.rm=T),0) #extract the summary data\n',sep='',file=zz)
					} else { 
						cat('        tt = tdata[,coi.in:(coi.in-',ttime,')] #extract a subset of the data\n',sep='',file=zz)
						cat('        for (ii in 1:12) {tt[,ii] = rowMeans(tt[,seq(ii,ncol(tt),12)],na.rm=TRUE)}; tt = tt[,1:12] #summarize to get averages\n',sep='',file=zz)
						cat('        tt = round(rowSums(tt,na.rm=T),0) #extract the summary data\n',sep='',file=zz)
					}
				} else if (voi=='bc15') {
					if (ttime<12) {
						cat('        tt = round((apply(tdata[,coi.in:(coi.in-',ttime,')],1,new.sd)/(rowMeans(tdata[,coi.in:(coi.in-',ttime,')],na.rm=T)+1))*100,0) #extract the summary data\n',sep='',file=zz)
					} else { 
						cat('        tt = tdata[,coi.in:(coi.in-',ttime,')] #extract a subset of the data\n',sep='',file=zz)
						cat('        for (ii in 1:12) {tt[,ii] = rowMeans(tt[,seq(ii,ncol(tt),12)],na.rm=TRUE)}; tt = tt[,1:12] #summarize to get averages\n',sep='',file=zz)
						cat('        tt = round((apply(tt,1,new.sd)/(rowMeans(tt,na.rm=T)+1))*100,0)  #extract the summary data\n',sep='',file=zz)
					}	
				} else if (voi=='bc16'){
					cat('        tt = tdata[,coi.in:(coi.in-',ttime,')] #extract a subset of the data\n',sep='',file=zz)
					if (ttime>12) cat('        for (ii in 1:12) {tt[,ii] = rowMeans(tt[,seq(ii,ncol(tt),12)],na.rm=TRUE)}; tt = tt[,1:12] #summarize to get averages\n',sep='',file=zz)
					cat('        tvals=c(1:12,1:12) #just some temp object to extract circular values (e.g. months 11,12,1 for 3 month FI)\n',sep='',file=zz)
					cat('        prec.monthly.mean = t.data = tt #extract only the info for the year\n',sep='',file=zz)
					cat('        for(ii in 1:12) {t.data[,ii] = rowSums(prec.monthly.mean[,c(tvals[(ii+12):(ii+12-2)])],na.rm=T)} #cycle through and get the quarterly sums\n',sep='',file=zz)
					cat('        tt = round(apply(t.data,1,new.max),0) #extract the summary data\n',sep='',file=zz)
				} else if (voi=='bc17'){
					cat('        tt = tdata[,coi.in:(coi.in-',ttime,')] #extract a subset of the data\n',sep='',file=zz)
					if (ttime>12) cat('        for (ii in 1:12) {tt[,ii] = rowMeans(tt[,seq(ii,ncol(tt),12)],na.rm=TRUE)}; tt = tt[,1:12] #summarize to get averages\n',sep='',file=zz)
					cat('        tvals=c(1:12,1:12) #just some temp object to extract circular values (e.g. months 11,12,1 for 3 month FI)\n',sep='',file=zz)
					cat('        prec.monthly.mean = t.data = tt #extract only the info for the year\n',sep='',file=zz)
					cat('        for(ii in 1:12) {t.data[,ii] = rowSums(prec.monthly.mean[,c(tvals[(ii+12):(ii+12-2)])],na.rm=T)} #cycle through and get the quarterly sums\n',sep='',file=zz)
					cat('        tt = round(apply(t.data,1,new.min),0) #extract the summary data\n',sep='',file=zz)
				}
				cat('        tasc = base.asc; tasc[cbind(pos$row,pos$col)] = tt #replace the base values for the data\n',sep='',file=zz)
				cat('        write.asc(tasc,paste(yy,mm,"/',voi,'.',ttime+1,'m.asc",sep="")) #write out the  ascii grid file\n',sep='',file=zz)
				cat('    } ; cat("\\n")\n',sep='',file=zz)
				cat('} ; cat("\\n")\n',sep='',file=zz)

			close(zz)
			#submit the job
			tt = ifelse(UQ,paste('qsub -l select=1:ncpus=8:NodeType=medium -l walltime=100:00:00 -A q1086 ',voi,ttime,'.sh',sep=''),paste('qsub -l nodes=1:ppn=8 ',voi,ttime,'.sh',sep=''))
			system(tt)
		}
	}
	system('sleep 60')
}


###########################################################################################
###########################################################################################
#create the csv files for the ascii grid files

#get a list of the directories that are going to be summarized
asc.dirs = list.files(work.dir)
#define the position file
pos = paste(base.data.dir,'summaries/Oz/base.positions.csv',sep='')

#go to the temporary pbs directory
pbs.dir = paste(base.dir,'tmp.pbs/',sep=''); system(paste('rm -rf ',pbs.dir,'*',sep='')); setwd(pbs.dir)

for (asc.dir in asc.dirs) {
	##create the sh file
	zz = file(paste(asc.dir,'.sh',sep=''),'w')
		cat('##################################\n',file=zz)
		cat('#!/bin/bash\n',file=zz)
		cat('source /etc/profile \n',file=zz)
		cat('module load R \n',file=zz) #load the necessary module				
		cat('cd $PBS_O_WORKDIR\n',file=zz)
		cat('R CMD BATCH --no-save ',asc.dir,'.R\n',sep='',file=zz)
		cat('##################################\n',file=zz)
	close(zz)
	#create the R summarizaing script
	zz = file(paste(asc.dir,'.R',sep=''),'w')
		cat('library(SDMTools)\n',file=zz)
		cat('#define directories\n',sep='',file=zz)
		cat('work.dir = "',work.dir,asc.dir,'/"; setwd(work.dir)\n',sep='',file=zz)
		cat('#read in necessary data\n',file=zz)
		cat('pos = read.csv("',pos,'")\n',sep='',file=zz)
		cat('for (asc.file in list.files(,pattern="\\\\.asc")) {\n',sep='',file=zz)
		cat('    tasc = read.asc(asc.file)\n',sep='',file=zz)
		cat('    pos[gsub("\\\\.asc","",asc.file)] = tasc[cbind(pos$row,pos$col)]\n',sep='',file=zz)
		cat('}\n',sep='',file=zz)
		cat('write.csv(pos,"data.csv",row.names=FALSE)\n',sep='',file=zz)
	close(zz)
	#submit the job
	tt = ifelse(UQ,paste('qsub -l select=1:ncpus=1:NodeType=medium -l walltime=100:00:00 -A q1086 ',asc.dir,'.sh',sep=''),paste('qsub ',asc.dir,'.sh',sep=''))
	system(tt)
}

###########################################################################################
###########################################################################################
#convert asc files to mxe

#get a list of the directories that are going to be summarized
asc.dirs = list.files(work.dir)

#define files/folders
mxe.dir = gsub('monthly.asc','monthly.mxe',work.dir); dir.create(mxe.dir)
maxent = paste(base.data.dir,'summaries/maxent.jar',sep='')

#define the tmp pbs directory and clean it
pbs.dir = paste(base.dir,'tmp.pbs/',sep=''); system(paste('rm -rf ',pbs.dir,'*',sep='')); setwd(pbs.dir)

for (asc.dir in asc.dirs) {
	#create the out.dir
	out.dir = paste(mxe.dir,asc.dir,'/',sep=''); dir.create(out.dir,recursive=T)
	##create the sh file
	zz = file(paste(asc.dir,'.sh',sep=''),'w')
		cat('##################################\n',file=zz)
		cat('#!/bin/bash\n',file=zz)
		cat('source /etc/profile \n',file=zz)
		cat('cd $PBS_O_WORKDIR\n',file=zz)
		cat('java -mx1024m -cp ',maxent,' density.Convert ',work.dir,asc.dir,'/ asc ',out.dir,' mxe\n',sep='',file=zz)
		cat('##################################\n',file=zz)
	close(zz)
	#submit the job
	tt = ifelse(UQ,paste('qsub -l select=1:ncpus=1:NodeType=medium -l walltime=100:00:00 -A q1086 ',asc.dir,'.sh',sep=''),paste('qsub -m n -l nodes=1:ppn=1 ',asc.dir,'.sh',sep=''))
	system(tt); system('sleep 1')
}
